
module calc_eph
implicit none
! arrays for 11 objects
    ! 1- Mercury, 2- Venus  , 3- Earth+Moon,  4- Mars,  5- Jupiter 6- Saturn,
    ! 7- Uranus , 8- Neptune, 9- Pluto ,     10- Moon, 11- Sun.
    integer jdb(1:11)              !initial Julian date (int part)
    double precision djb(1:11)  !initial Julian date (frac part)
    double precision delta(1:11)!approximation subinterval
    integer ncoef(1:11)            !number of (polunomial degrees + free term)
    integer ndim(1:11)             !number of dimensions
    integer nbl(1:11)              !number of subintervals
    integer shift(1:11)            !shift in the coeef array
    double precision coef(1:Z'4000000')
                                   !array of coefficiens for 11 objects
    character(80) names(1:11)      !file names
    double precision, parameter:: rho = 82.3005676536174207 !const Earth/Moon + 1

contains
logical function InitTxt(dir)
!Read 11 ephemeride text files into arrays
!Input: path to the directory '\TEXT\'
!Output: logical
    character(*):: dir
    integer:: dir_l
    integer:: i, ios, coef_size, shift_, lfn = 77

    InitTxt = .TRUE.
    dir_l = len_trim(dir)
    do i = 1, 11
      names(i)(1:dir_l) = dir
    end do
    write(names(1)(dir_l + 1:80),'(A)') 'text/mercury.08t'
    write(names(2)(dir_l + 1:80),'(A)') 'text/venus.08t'
    write(names(3)(dir_l + 1:80),'(A)') 'text/earth_m.08t'
    write(names(4)(dir_l + 1:80),'(A)') 'text/mars.08t'
    write(names(5)(dir_l + 1:80),'(A)') 'text/jupiter.08t'
    write(names(6)(dir_l + 1:80),'(A)') 'text/saturn.08t'
    write(names(7)(dir_l + 1:80),'(A)') 'text/uranus.08t'
    write(names(8)(dir_l + 1:80),'(A)') 'text/neptune.08t'
    write(names(9)(dir_l + 1:80),'(A)') 'text/pluto.08t'
    write(names(10)(dir_l + 1:80),'(A)') 'text/moon.08t'
    write(names(11)(dir_l + 1:80),'(A)') 'text/sun.08t'

    shift_ = 1
    do i = 1, 11
        write(*,'(A60)') names(i)
        open(unit = lfn, file = names(i), status = 'old', iostat = ios)
        if(ios.NE.0) then
          write(*,*) "Not found"
          InitTxt = .FALSE.
          return
        end if
        write(*, *) '.....'
        read(unit = lfn, fmt='(I25)') jdb(i)
        read(unit = lfn, fmt='(F25.17)') djb(i)
        read(unit = lfn, fmt='(F25.17)') delta(i)
        read(unit = lfn, fmt='(I25)') ncoef(i)
        read(unit = lfn, fmt='(I25)') ndim(i)
        read(unit = lfn, fmt='(I25)') nbl(i)

        coef_size = ncoef(i) * ndim(i) * nbl(i)
        read(lfn, fmt='(F25.17)') coef(shift_: shift_ + coef_size - 1)
        shift(i) = shift_
        shift_ = shift_ + coef_size
        close(lfn)
        write(*,*) "Ok"
     end do
end function InitTxt

logical function InitBin(dir)
!Read 11 ephemeride binary files into arrays
!Input: path to the directory '\BIN\'
!Output: logical
!temporary variables to set up a correspondence between 26b title and global arrays
    character(*):: dir
    integer:: dir_l

    integer jdb_
    double precision    djb_
    double precision    delta_
    integer ncoef_
    integer ndim_
    integer nbl_
    integer shift_

    integer title(8)  !first 32 bytes of the file
    equivalence (title(1),jdb_)
    equivalence (title(2),djb_)
    equivalence (title(4),delta_)
    equivalence (title(6),ncoef_)
    equivalence (title(7),ndim_)
    equivalence (title(8),nbl_)

    integer:: i, j, ios, coef_size, lfn = 77

    InitBin = .TRUE.
    dir_l = len_trim(dir)
    do i = 1, 11
      names(i)(1:dir_l) = dir
    end do
    write(names(1)(dir_l + 1:80),'(A)') 'BIN/mercury.08b'
    write(names(2)(dir_l + 1:80),'(A)') 'BIN/venus.08b'
    write(names(3)(dir_l + 1:80),'(A)') 'BIN/earth_m.08b'
    write(names(4)(dir_l + 1:80),'(A)') 'BIN/mars.08b'
    write(names(5)(dir_l + 1:80),'(A)') 'BIN/jupiter.08b'
    write(names(6)(dir_l + 1:80),'(A)') 'BIN/saturn.08b'
    write(names(7)(dir_l + 1:80),'(A)') 'BIN/uranus.08b'
    write(names(8)(dir_l + 1:80),'(A)') 'BIN/neptune.08b'
    write(names(9)(dir_l + 1:80),'(A)') 'BIN/pluto.08b'
    write(names(10)(dir_l + 1:80),'(A)') 'BIN/moon.08b'
    write(names(11)(dir_l + 1:80),'(A)') 'BIN/sun.08b'

    title = 0
    shift_ = 1
    do i = 1, 11
        write(*,'(A60)') names(i)
        open(unit = lfn, file = names(i), form = 'unformatted', &
        &    access = 'direct', recl = 4, status = 'old', iostat = ios)
        if(ios.NE.0) then
          write(*,*) "Not found"
          InitBin = .FALSE.
          return
        end if
        write(*, *) '.....'
        do j = 1, 8
          read(unit = lfn, rec = j) title(j)
        end do
        close(lfn)
        open(unit = lfn, file = names(i), form = 'unformatted', &
        &    access = 'direct', recl = 8, status = 'old', iostat = ios)
        jdb(i) = jdb_
        djb(i) = djb_
        delta(i) = delta_
        ncoef(i) = ncoef_
        ndim(i) = ndim_
        nbl(i) = nbl_

        coef_size = ncoef_ * ndim_ * nbl_
        shift(i) = shift_
        do j = 1, coef_size
          read(unit = lfn, rec = 4 + j) coef(j + shift_ - 1)
        end do
        shift_ = shift_ + coef_size
        close(lfn)
        write(*,*) "Ok"
     end do
end function InitBin

subroutine calc_EPM(pl_num, centr_num, jd, dj, xv, vv, status)
!Calculates coordinates and velocities of the oblect 'PL_NUM'
!with respect to the point 'CENTR_NUM'
!The numbering convention for 'PL_NUM' and 'CENTR_NUM':
!
!                1 = MERCURY           8 = NEPTUNE
!                2 = VENUS             9 = PLUTO
!                3 = EARTH            10 = MOON
!                4 = MARS             11 = SUN
!                5 = JUPITER          12 = SOLAR-SYSTEM BARYCENTER
!                6 = SATURN           13 = EARTH-MOON BARYCENTER
!                7 = URANUS
    !Input
    integer:: pl_num    !number of object
    integer:: centr_num !number of center
    integer:: jd        !Julian date (integral part) at which interpolation is wanted
    double precision:: dj           !Julian date (frac part)
    !Output
    double precision:: xv(1:3)      !object positions
    double precision:: vv(1:3)      !object velocities
    logical:: status       !status of output data

    !Internal
    double precision:: xv_bar(1:3)      !barycentric positions of the point 'CENTR_NUM'
    double precision:: vv_bar(1:3)      ! -"- velocities  -"-

    status = .TRUE.
    xv = 0
    vv = 0

    if((pl_num.EQ.centr_num).OR.(pl_num.EQ.12)) then
      return
    end if

    !Moon(geo)
    if((centr_num.EQ.3).AND.(pl_num.EQ.10)) then
       call calc(pl_num, jd, dj, xv, vv, status)
       return
    end if

    !barycentric
    if(centr_num.EQ.12) then
       call calc_bar(pl_num, jd, dj, xv, vv, status)
       return
    end if

    !general case
    call calc_bar(pl_num, jd, dj, xv, vv, status)
    if(.NOT.status) return
    call calc_bar(centr_num, jd, dj, xv_bar, vv_bar, status)
    if(.NOT.status) return
    xv = xv - xv_bar
    vv = vv - vv_bar

end subroutine calc_EPM

subroutine calc_bar(pl_num, jd, dj, xv, vv, status)
!Calculates barycentric coordinates and velocities of the following objects
! 1- Mercury, 2- Venus  , 3- Earth,  4- Mars,  5- Jupiter
! 6- Saturn,  7- Uranus , 8- Neptune, 9 - Pluto,
! 10- Moon, 11- Sun, 13 - Earth + Moon
    !Input
    integer:: pl_num    !number of object
    integer:: jd        !Julian date (integral part) at which interpolation is wanted
    double precision:: dj           !Julian date (frac part)
    !Output
    double precision:: xv(1:3)      !object positions
    double precision:: vv(1:3)      !object velocities
    logical:: status       !status of output data
    !Internal
    double precision:: xv_add(1:3)      !additional positions
    double precision:: vv_add(1:3)      !additional velocities

    if(pl_num.EQ.3) then !Earth
        !Earth(bar) = (Earth + Moon)(bar) - Moon(geo)/rho
        call calc(3, jd, dj, xv, vv, status)
        if(.NOT.status) return
        call calc(10, jd, dj, xv_add, vv_add, status)
        if(.NOT.status) return
        xv_add = xv_add / rho
        vv_add = vv_add / rho
        xv = xv - xv_add
        vv = vv - vv_add
        return
    end if

    if(pl_num.EQ.10) then !Moon
        !Moon(bar) = Earth(bar) + Moon(geo)
        call calc(3, jd, dj, xv, vv, status)
        if(.NOT.status) return
        call calc(10, jd, dj, xv_add, vv_add, status)
        if(.NOT.status) return
        xv_add = xv_add / rho
        vv_add = vv_add / rho
        xv = xv - xv_add
        vv = vv - vv_add
        !Moon(geo)
        call calc(10, jd, dj, xv_add, vv_add, status)
        if(.NOT.status) return
        xv = xv + xv_add
        vv = vv + vv_add
        return
    end if

    if(pl_num.EQ.13) then !Earth + Moon
        call calc(3, jd, dj, xv, vv, status)
        return
    end if

    !neither Earth, nor Moon
    call calc(pl_num, jd, dj, xv, vv, status)
end subroutine calc_bar

subroutine calc(pl_num, jd, dj, xv, vv, status)
!Calculates coordinates and velocities of the following oblects:
! 1- Mercury(bar), 2- Venus(bar)  , 3- Earth+Moon(bar),  4- Mars(bar),  5- Jupiter(bar)
! 6- Saturn(bar),  7- Uranus(bar) , 8- Neptune(bar), 9 - Pluto(bar) ,
! 10- Moon(geo), 11- Sun(bar).
    !Input
    integer:: pl_num    !number of object
    integer:: jd        !Julian date (integral part) at which interpolation is wanted
    double precision:: dj           !Julian date (frac part)
    !Output
    double precision:: xv(1:3)      !object positions
    double precision:: vv(1:3)      !object velocities
    logical:: status       !status of output data

    !local variables
    integer:: jdb_          ! initial Julian date (int part)
    double precision::    djb_          ! initial Julian date (frac part)
    double precision::    delta_        ! approximation subinterval
    integer:: nbl_          ! number of subintervals
    integer:: ncoef_        ! number of (polunomial degrees + free term)
    integer:: ncoefm        ! number of (polunomial degrees) ncoef_ - 1
    integer:: ndim_         ! number of dimensions
    integer:: shift_        ! position of coeffiients of current object in the array coef

    double precision:: td               ! Julian in [initial date, final date] -> td in [-1, 1]
    integer:: cur_bl        ! approximation interval for current Julian date
    double precision:: tt(32), xx(32)   ! arrays of values and integrals of Chebyshev polynomials
    integer:: coef_ptr         ! position of the current interval in the array coef
    integer:: i                ! counter

    jdb_ = jdb(pl_num)
    djb_ = djb(pl_num)
    nbl_ = nbl(pl_num)
    ndim_ = ndim(pl_num)
    delta_ = delta(pl_num)
    ncoef_ = ncoef(pl_num)
    shift_ = shift(pl_num)

    ncoefm = int(ncoef_ - 1, 1)

    td = ((jd - jdb_) + (dj - djb_)) / delta_
    cur_bl = int(td)
    td = td - cur_bl
    td = 2.0 * td - 1.0
    cur_bl = cur_bl + 1
    if((cur_bl.LE.0).OR.(cur_bl.GT.nbl_)) then
      status = .FALSE.
      return
    end if
    shift_ = shift_ + (cur_bl - 1) * (ncoef_ * ndim_)
    call calc_cheb_pol(ncoefm, td, tt, xx)
    do i = 1, ndim_
       coef_ptr = shift_ + ncoef_ * (i - 1)
       call calc_pos(ncoefm, coef_ptr, xx, delta_, xv(i))
       call calc_vel(ncoefm, coef_ptr, tt, vv(i))
    end do
end subroutine calc

subroutine calc_cheb_pol(ncoefm, td, tt, xx)
! Calculates values of Chebyshev's polynomials along with their
! integrals.
       ! Input
       integer::ncoefm            ! polynomial power
       double precision:: td                  ! argument [-1 .. 1]
       !Output
       double precision:: tt(:)               ! array of values of Chebyshev polynomials
       double precision:: xx(:)               ! array of integrals

       double precision:: d                   ! 2.0 * td
       integer:: i, j, k          ! counters
       logical:: flag                ! (-1)^n

       tt(1) = 1.0
       tt(2) = td
       d = 2.0 * td
       do i = 3, ncoefm + 1
           tt(i) = tt(i - 1) * d - tt(i - 2)
       end do
       xx(1) = td
       xx(2) = (tt(3) + tt(1)) * 0.25
       do i = 3, ncoefm
           xx(i) = 0.5 * (tt(i + 1) / i - tt(i - 1) / (i - 2))
       end do
       j = 0
       flag = .false.
       i = 4
       do  k = 1, ncoefm / 2 - 1
           j = j + 1
           d = 0.25 / j + 0.25 / (j + 1)
           flag = .not.flag
           if (flag) d = -d
           xx(i) = xx(i) + d
           i = i + 2
       end do
end subroutine calc_cheb_pol

subroutine calc_vel(ncoefm, c_ptr, tt, v)
! Calculates value of velocity's component.
       !Input
       integer:: ncoefm   ! number of coefficients
       integer:: c_ptr    ! position of the current interval in the array coef
       double precision::  tt(:)      ! array of values of Chebyshev's polynomials
       !Output
       double precision:: v           !velocity's component

       integer:: i
       v = 0.0
       do i = ncoefm, 2, -1
         v = v + coef(c_ptr + i - 1) * tt(i)
       end do
       v = v + coef(c_ptr)
end subroutine

subroutine calc_pos(ncoefm, c_ptr, xx, dlt, x)
!Calculates value of object coordinate
       !Input
       integer:: ncoefm     ! number of coefficients
       integer:: c_ptr      ! position of the current interval in the array coef
       double precision:: xx(:)         ! array of integrals of Chebyshev's polynomials
       double precision:: dlt           ! approximation interval
       !Output
       double precision:: x             ! object coordinate

       integer:: i, j
       double precision:: x0            ! free term, c(c_ptr + ncoefm)

       x = 0.0
       do j = 1, ncoefm
         i = ncoefm - j + 1
         x = x + coef(c_ptr + i - 1) * xx (i)
       end do
       x0 = coef(c_ptr + ncoefm)
       x = 0.5 * dlt * x + x0
end subroutine calc_pos

! ===============================================================

subroutine InitEpm(bTxt, dir)
    logical     bTxt
    logical     status
    character(*):: dir
    !--------------------------------
    if(bTxt) then
        ! To read text ephemeride files call function  InitTxt(dir)
        ! if 'DIR' is empty string search in the work directory\Text\
        status = InitTxt(dir)
    else
        ! To read binary ephemeride files call function  InitBin(dir)
        ! if 'DIR' is empty string search in the work directory\BIN\
        status = InitBin(dir)
    end if

    if(.NOT. status) then
      write(*, *) 'error loading EPM ephemeride'
      return
    end if
end subroutine InitEpm

Subroutine InitEphemTimeRangeEPM()
    implicit none
    REAL*8 EthenTimeRange(3)
    COMMON /EthenTimeRange/EthenTimeRange

      integer:: jdb_          ! initial Julian date (int part)
    double precision::    djb_          ! initial Julian date (frac part)
    double precision::    delta_        ! approximation subinterval
    integer:: nbl_          ! number of subintervals
    integer:: ncoef_        ! number of (polunomial degrees + free term)
    integer:: ncoefm        ! number of (polunomial degrees) ncoef_ - 1
    integer:: ndim_         ! number of dimensions
    integer:: shift_        ! position of coeffiients of current object in the array coef
    !integer:: pl_num
    integer:: pl_num
    !--------------------------------

    pl_num = 3 ! assume Earth
    jdb_ = jdb(pl_num)
    djb_ = djb(pl_num)
    nbl_ = nbl(pl_num)
    ndim_ = ndim(pl_num)
    delta_ = delta(pl_num)
    ncoef_ = ncoef(pl_num)
    shift_ = shift(pl_num)

    ncoefm = int(ncoef_ - 1, 1)

    EthenTimeRange(1) = jdb_ + djb_
    EthenTimeRange(2) = EthenTimeRange(1) + delta_ * nbl_
    EthenTimeRange(3) = delta_;
end Subroutine InitEphemTimeRangeEPM


subroutine InitGMsEPM()
    implicit none

    real*8 GM(11)
    COMMON /GM/GM

    real*8 kmAE
    common /kmAE/kmAE
!--------------------------------

    kmAE = 149597870.696547

    GM(1) =     49125.4957D-15    ! Mercury
    GM(2) =    724345.2333D-15  ! Venus
    GM(3) =    899701.1347D-15  ! Earth + Moon
    GM(4) =     95495.4869D-15  ! Mars
    GM(5) = 282534584.0787D-15  ! Jupiter
    GM(6) =  84597060.7310D-15  ! Saturn
    GM(7) =  12920248.2576D-15  !    Uranus
    GM(8) =  15243591.0921D-15  ! Neptune
    GM(9) =      2166.8443D-15  ! Pluto
    GM(10)=         0.0D0         ! Moon
    GM(11)= 0.01720209895D0 * 0.01720209895D0  ! Sun

    GM(10)= GM(3) / rho     ! Moo
    GM(3) = GM(3) - GM(10)            ! Earth
end subroutine InitGMsEPM

end module calc_eph
